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Abstract —One of the goals of computer-aided surgery is 
to match intraoperative data to preoperative images of the 
anatomy and add complementary information that can facilitate 
the task of surgical navigation. In this context, mechanical 
palpation can reveal critical anatomical features such as arteries 
and cancerous lumps which are stiffer that the surrounding tis¬ 
sue. This work uses position and force measurements obtained 
during mechanical palpation for registration and stiffness 
mapping. Prior approaches, including our own, exhaustively 
palpated the entire organ to achieve this goal. To overcome the 
costly palpation of the entire organ, a Bayesian optimization 
framework is introduced to guide the end effector to palpate 
stiff regions while simultaneously updating the registration 
of the end effector to an a priori geometric model of the 
organ, hence enabling the fusion of intraoperative data into the 
a priori model obtained through imaging. This new framework 
uses Gaussian processes to model the stiffness distribution 
and Bayesian optimization to direct where to sample next 
for maximum information gain. The proposed method was 
evaluated with experimental data obtained using a Cartesian 
robot interacting with a silicone organ model and an ex vivo 
porcine liver. 

1. INTRODUCTION 

Surgeons performing minimally invasive surgery (MIS) 
offer their patients a shorter recovery time and reduced 
pain at a cost of increased technical difficulty. One of the 
drawbacks of MIS is the loss of direct sensory feedback. 
This loss can impede the detection and use of surface 
and stiffness features which can help the surgeon find 
correspondence between the intraoperative scene and the 
preoperative imaging data. Computer-aided surgery (CAS) 
was introduced to provide surgeons performing MIS 
with functional and geometric information that can aid 
intraoperative navigation and execution of the preoperative 
plans. Mechanical exploration by palpation and tissue 
manipulation can provide complementary information about 
geometric constraints [1], [2], tissue characteristics [3] 
and variation in stiffness throughout the organ [4], thus 
augmenting the information obtained through conventional 
image-based CAS. In this paper, we introduce tools for 
efficient information-guided exploration of an organ using 
mechanical palpation and integration of information, such 
as stiffness, to an a priori model through registration. 
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Mechanical palpation can facilitate the localization of 
arteries and other anatomically critical features which 
visually cannot be detected [5]-[7]. Several groups 
previously proposed heuristic algorithms for autonomous 
exploration [8] and segmentation of stiff features in 
compliant environments using classifiers [9]. These 
methods rely on the robot visiting all locations on 
a discretized grid on the organ’s surface. Grid maps 
comes with the assumption of independence between 
the grid points ignoring the structural dependency in the 
environment. Accurate registration of the surgical tool to 
the coordinate frame of the preoperative model is one 
of the primary goals in incorporating intraoperative data, 
such as stiffness, into an a priori geometric model of the 
anatomy. Aforementioned previous works on palpation 
and exploration are not concerned with incorporating the 
stiffness information into the preoperative geometric data, 
a task that can help the surgeon in correlating preoperative 
models, like CT scans, with the current intraoperative scene. 

To avoid the costly palpation of the entire organ and 
enable fusion of the stiffness information into an a priori 
geometric model of the anatomy, we introduce a method 
based on Bayesian optimization that minimizes the amount 
of probing required to reveal stiff features and registers 
the tool to the a priori model. The proposed algorithm 
is advantageous because it only visits regions that bring 
information gain, contrary to searching a discretized grid on 
the entire surface of the organ to find stiff features. 

In the Bayesian formulation, a Gaussian process (GP) [10] 
is used to define a prior over the unknown stiffness 
distribution. Gaussian Processes provide a probabilistic 
description of the stiffness map and captures the variance 
of the stiffness distribution which helps guide the probing 
towards unexplored regions. The stiff regions correspond 
to regions near local maxima of the stiffness distribution 
and the Bayesian optimization finds the maxima of this 
unknown stiffness distribution by directing the probing to 
points that would result in maximum information gain in 
predicting the stiff regions. In a complementary effort, our 
collaborators at Johns Hopkins University are exploring 
different GP formulations to concurrently estimate surface 
geometry and stiffness for model reconstruction, using 
continuous palpation motion. 

We first introduce GP and Bayesian optimization in Sec¬ 
tion |n| followed by the description of simultaneous registra- 


tion and stiffness estimation. The experimental evaluations 
and results are given in Section m 


II. BACKGROUND 


A. Gaussian Processes 


A stochastic process is a collection of random variables, 
{Y : X G A'}, indexed by elements from a set A’, known as 
the index set. A Gaussian process is a stochastic process 
such that any finite subcollection of random variables has 
a multivariate Gaussian distribution [11]. A GP , / ^ 
k), is fully specified by its mean function /i : A' —> M 
and a covariance function k : A! x A! ^ M+. 

Intuitively, we can think of GP as a distribution over func¬ 
tions. Each random variable, Yi, in a GP’s collection is the 
distribution of function values at a point xi ^ X. Gaussian 
processes can be used for regression and to make predictions 
at a new point x* G A by defining a prior over functions. 
Given a set of n observed inputs x = [xi, X 2 ,..., 
and corresponding outputs Y = [Yi, I 2 , • • •, the 

random variables Y are Gaussian distributed with mean 
/i(x 2 ),..., and covariance matrix K whose 

elements are defined by a covariance function, k{xi,Xj) 
where i, j G [1, ...n], that defines the covariance between Yi 
and Yj . A commonly used covariance function is the squared 
exponential kernel defined as 

k{xi, Xj) = Gf exp (^— 

where crj is the variance of the process used as a scaling 
factor and £ is the length-scale of the kernel. 

If we want to make predictions at a new set of m points 
X*, the joint distribution is given as 
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In machine learning literature, x is called the training set and 
cc* is called the prediction or test set [10]. For simplicity, the 
prior mean function is generally assumed to be zero mean. 
The conditional (predictive) distribution Y* can be computed 
using the conditioning rule of multivariate Gaussian distri¬ 
butions 


p{Y,\Y) r. M{K,K-^Y,K,, - K,K-^Kl). (3) 

where is the m x n training-prediction set covariance 
X(x*,x) and is the m x m prediction set covariance 
matrix x*). 

We employ GP to model the stiffness distribution of the 
organ. The values of Y are the stiffness values associated 
with the probed points. The position of the points which are 
probed form the prediction set, x, while x* is the prediction 
grid which spans the surface of the organ. Note that GP 
is continuous and the prediction grid is only used to plot 
the stiffness distribution for visualization. By using GP, we 
assume the stiffness distribution changes smoothly across 
the organ and this smoothness is defined by the choice 
of the covariance function. In our formulation, we use a 


local deformation model for stiffness estimation. Therefore, 
a kernel with a fixed length-scale that is on the same order 
of the size of the palpation tool’s tip is an effective choice. 
We use the squared exponential kernel with df = 1 and 
i = 3mm. The length-scale, £, determines how close two 
points have to be for the observation at those points to be 
correlated. 


TABLE I 
Notation 


Symbol Description 

X Coordinates of the probed points in the training set 

£c* Coordinates of the grid points in the prediction set 

Y The output at probed points x 
V* The predicted output at cc* 

(yu, a) Mean and variance of the predictive distribution 
[•]^ Entity expressed in tool’s reference frame 
[•]^ Entity expressed in CAD model’s reference frame 
[•]o Initial value of the entity 
T Homogeneous transformation matrix 
n Normal vector 

p Coordinates of sensed position 

F Magnitude of sensed normal force 
c Stiffness at the probed point x 
cj) CAD model with triangle faces and vertices 


B. Exploration and Exploitation with Bayesian Optimization 

Bayesian optimization is a powerful framework for global 
optimization of black-box functions [12]. It is most beneficial 
when the function does not have a closed-form expression 
and obtaining observations from the function is expensive. 
Bayesian optimization allows for prior beliefs about an 
unknown function to be updated via a posterior. Stiffness 
distribution of an organ can be thought of the unknown 
function we want to optimize whose maxima correspond 
to the stiff features. In the Bayesian framework, we use 
GP to define a prior over the stiffness distribution. The 
sequential nature of the Bayesian optimization can help guide 
the sampling of the continuous search space. Sequential sam¬ 
pling requires selecting an acquisition function, also known 
as the utility function [13]. Acquisition functions use the 
mean, /i(cc), and variance, a{x),of the predictive distribution 
posterior to compute a function which shows the most likely 
locations of the global maximum.Acquisition functions such 
as probability of improvement [14], expectation improvement 
(El) [12], and upper-confidence based methods [15] have 
been developed to find the global maximum and to balance 
exploration with exploitation. El provides global exploration 
of the search space and local exploitation. The El acquisition 
function can be evaluated analytically and is given as [12]: 


EI{x) 


({n{x) — y+)<E(Z) + a{x)(l){Z) if a{x) > 0 
^0 if a{x) = 0 
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is the output at x~^ which is the current maximum of 
the sampled points. and ^(.) are the PDF and CDF of 
the standard normal distribution respectively. 

Fig. illustrates an expected improvement scenario for a 
ID example. Expectation improvement is generally used to 
find the global maximum of the unknown function and the 
search is terminated once a desired improvement is achieved. 
Note that after finding the global maximum, the algorithm 
still continues exploring other local maxima in the next 
iterations in an effort to reduce uncertainty and to find stiff 
features. In addition to using El, we do pure exploration 
after every 5 samples and select a random point that has 90% 
uncertainty. Such exploration is advantageous when there are 
multiple stiff features and the initial samples do not include 
points near stiff regions. Preoperative information such as the 
size of the tumor or the width of the artery is useful to decide 
on the density of the samples in the initial set, however this 
is not the explored in this work. Interested readers can refer 
to [16] for discussions on the effect of the initial sample 
size in GP predictions. 


n=4 




Fig. 1. A 1-D example that starts with an initial training set of 4 points. 
Red triangle shows the maximum of the El acquisition function which is 
the point that should be probed next. 
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Fig. 2. Two examples that show the advantage of using El as a sampling 
strategy: (a) and (b) show the ground truth stiffness map for a simulated data 
with an artery and simulated data with multiple stiff features, respectively; 
(c) and (d) show the maps obtained using uniformly sampling with 100 
points; (e) and (f) show the corresponding maps obtained using El with 
100 samples. 


Fig. demonstrates the stiffness map obtained using 
uniform sampling in comparison with the map sampled 
using El for two different cases. With El we can acquire 
useful information about the stiffness distribution compared 
to uniformly sampling the surface of the organ for the 
examples shown. 


C. Registration and Stiffness Estimation 

Our group has previously developed a method for simulta¬ 
neously estimating the registration and stiffness distribution 
over the surface of a fiexible environment using a Kalman 
filtering approach called CARE [17] and a more recent 
model update method, called Complementary Model Update 
(CMU), that decouples stiffness estimation from registration, 
resulting in a more robust implementation. Similar to CARE, 
the CMU uses the force and position information obtained 
by interaction of the surgical tool with the organ to estimate 
the local stiffness and to register the organ to its preoperative 
model. It is assumed that the local surface deformations are 
only due to physical interaction of the surgical tool with the 
organ. We also assume that the surface of the organ is smooth 
and frictionless, thus the applied force is along the surface 
normal and increases with depth. Registration is performed 
by finding the transformation that takes the probed points 
defined in the tool frame to the corresponding points on the 
preoperative model of the anatomy. The preoperative model 
is a computer-aided design (CAD) model in the form of a 
triangular mesh. 

The stiffness at a probed point is estimated using a best 
line fit between the relative sensed deformation depths and 
sensed forces: 


Ci=L (ll(pf)i - ip^)4 , {{Fp)i - {F^)i)) (6) 

where Ci is the stiffness of the probed point, {p^)i and 
{p^)i are the coordinates of two distinct sensed positions 
expressed in the tool reference frame, R, corresponding to the 
probed point on the surface of the organ and {F^)i 

are the corresponding magnitude of sensed forces. The 
registration estimate is obtained using the estimated stiffness, 
sensed position and magnitude of the sensed force: 


T = argmin 


i=l 


P 


P 
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where T G SE{3) is the homogeneous transformation that is 
to be estimated, pf is the location of the probed point in the 
model’s reference frame, given by C, and the surface normal 
at the probed point is denoted by nf. More information 
about the CMU is provided in Appendix. 


III. EXPERIMENTS AND RESULTS 
A. Experimental setup 

A Cartesian robot with an open architecture controller 
was used to evaluate the framework proposed in this paper 
as shown in Eig. |3(a)| The robot is equipped with an ATI 
Nano43 6-axis force sensor at the robot end-effector (EE) 
and is capable of executing the hybrid motion/force tasks 
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Fig. 3. (a) Experimental setup used for ex vivo organ experiments. The 

tool frame, R, is located at the end effector,(b) Contact location and surface 
norm estimation 


described in Khatib [18]. We assume that the robot end 
effector has been positioned above the organ. We do not 
assume any prior knowledge of the local surface normal. 

Experiments were conducted with a silicon phantom organ 
and an ex vivo porcine liver. In the silicon phantom exper¬ 
iment the top surface of the organ was lubricated and in 
the ex vivo experiment the organ surface was hydrated to 
reduce friction during probing. In both experiments, a target 
region was defined by the user. This region was then used 
as a reference to generate a uniformly distributed grid map 
with uniform spacing in the x-y plane of the tool’s reference 
frame, R. For a particular reference probing location, the 
following procedures were repeated automatically to obtain 
the force/deformation profile data. 

1) Making high force contact: The robot controller is given 
a desired probing location Xp and a force magnitude 
of 0.5N. The hybrid force/motion controller decouples 
the combined commands into compatible (orthogonal) 
force/motion commands that direct the robot to achieve 
a desired position in the x-y plane and a stable contact 
force with the environment along the z axis. 

2) Estimating surface norm: The contact location and 
surface normal estimation is shown in Fig. |3(b)| The 
surface normal, n, is computed using the force sensed 
from the environment, n = fs/||fs||, assuming the 
surface friction is negligible. 

3) Finding low force surface contact point: In this step, 
the robot first retreats swiftly away from the surface 
and then moves towards the surface, along the surface 
normal this time, to find the zero force intersection. 
An offset is applied from the robot EE to obtain the 
estimated contact point as Xcont = ^ee — where r is 
the radius of the robot end effector ball. A 9mm radius 
probe was used in the silicon examples and a 6.3mm 
radius probe was used in the ex vivo organ experiment. 

4) Probing and recording: The robot is commanded in po¬ 
sition control mode along the estimated surface normal 
direction with 0.3mm increments until 3mm probing 
depth. Hence, there are 10 position measurements, p, 
and 10 force measurements, F, for each point we probe. 


B. Bayesian Optimization Guided Probing 

To evaluate Bayesian optimization guided probing, we 
simulated experiments using the experimentally collected 


data. A block diagram description of the probing method is 
shown in Fig.|^ Prediction of the stiffness distribution is car¬ 
ried out in the tool frame, R. Initially, the actual registration 
is unknown, hence we do not know where the probed points 
correspond to in the preoperative CAD model. We assume 
palpation is carried out inside a region of interest (ROI). It is 
assumed that the initial set of samples, Xq^, include points 
on the boundary of the ROI as well as uniformly distributed 
points inside the ROI. The training set consists of previously 
probed points, where i = l,2,...,n. We use gridfit 
function [19] to interpolate the previously probed points to 
form a dense grid inside the ROI, where j = 1,2, ...m, 
to make predictions using GP. This grid is used to estimate 
the stiffness distribution for visualization. The stiffness value 
at a probed point, q, is estimated, in our case by the 
CMU, and corresponds to at Xi^ in the GP formulation. 
Based on the posterior of the predictive distribution given 
by /i and a the point at which the expectation improvement 
takes a maximum value is selected to be the next palpation 
point, As a new point is palpated, it is added to 

the training set and the prediction grid is regenerated. We 
use the updated registration estimate, T, to transform the 
probed points and their associated stiffness values to the 
corresponding points on the CAD model. This procedure 
enables displaying experimentally collected stiffness data on 
the preoperative CAD model. 



Fig. 4. Block diagram description of the Bayesian optimization guided 
probing 


C. Evaluation 

The proposed method was evaluated with various stiffness 
distributions and in the presence of measurement noise. We 
test the algorithm for four different scenarios: 

1. Simulated data of an organ with a multimodal stiffness 
distribution. 

2. Simulated data of an organ with a perturbed multimodal 
stiffness distribution and sensor noise. 

3. Silicon phantom organ with an embedded mock artery. 

4. Ex vivo porcine liver with a stiff inclusion. 

The goal of Example 1 is to demonstrate that expectation 
improvement is effective in finding all the local maxima and 
not just the global maximum. We start with an initial set 
of 19 samples and terminate the palpation after 100 points. 
Fig. I^a) shows where we think the position of all the probed 
points are based on the initial registration guess and their 
registered position estimated by the CMU algorithm. The 





































registered position of the probed points (deformed points) lie 
below the surface of the CAD model as expected; validating 
that the registration estimate is accurate. Fig. [^b) and (c) 
show the ground truth stiffness map and the predicted stiff¬ 
ness map, respectively. The predicted stiffness map captures 
the stiff features present in the ground truth stiffness map. 






Fig. 5. The algorithm starts with an initial set of of 19 points: 4 corners 
and 15 uniformly spaced points. The results are shown for 119 probed 
points. Example 1: (a) Registration results, (b) Ground truth stiffness map, 
(c) Estimated stiffness map 

Example 2 shows the effect of noisy sensor measurements 
in the prediction of the stiffness distribution. The ground 
truth stiffness map of Example 2 was obtained by perturbing 
the stiffness distribution for Example 1. An artificial sensor 
noise with (/ix^cTx) = (0, 0.3mm) was added to the sensed 
position and (/iF,cr^) = (0, O.IN) was added to the sensed 
force to simulate a more realistic scenario. Fig. [^a) shows 
the registration results. Fig. [^b) and (c) show the ground 
truth stiffness map and the predicted stiffness map, respec¬ 
tively. In the presence of noisy sensor measurements, the 
algorithm still reveals the stiff features and the registration 
parameters converge to the correct values. 

The silicon organ used in Example 3 and the ex vivo 
organ used in Example 4 are shown in Fig. [Tfa) and (b), 
correspondingly. The ground truth stiffness distribution for 
Example 3 and for the ROI in Example 4 were obtained 
by interpolating the experimental data at the grid locations 
shown in Fig. [tJc) and (d) and are shown in Fig. |7je) and 
(f), respectively. We emphasize that the organ was discretized 
only to palpate each grid point with the robot for the purpose 
of generating a ground truth stiffness map of the organ 
to test our algorithm. Fig. [^b) shows that sensed force is 
proportional to the depth of probing for the three locations 
shown on the CAD model of the liver in Fig. [^a), validating 
that the linear stiffness is a valid assumption for 3mm 
probing depth. 

Fig. shows the successful registration and the estimated 
stiffness map that reveals the mock artery in the silicon 
organ. Fig. shows the registration result and the position 


Fig. 6. Example 2: (a) Registration results, (b) Ground truth stiffness map, 
(c) Estimated stiffness map. The algorithm starts with an initial set of 19 
points: 4 corners and 15 uniformly spaced points. The results are shown for 
119 probed points. 
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Fig. 7. (a) Silicon phantom organ (b) ex vivo porcine liver with an inclusion 
sutured inside the organ, (c) 619 points were probed on the organ to generate 
a ground truth for the silicon organ, (d) 196 points were probed on the 
ex vivo organ to generate a ground truth stiffness map, (e) Stiffness map 
used as the ground truth for the silicon organ, (f) Stiffness map used as the 
ground truth for the ex vivo organ 



Fig. 8. (a) CAD model of the organ showing three locations with a square, 
star and a circle, (b) Force vs probing depth for the three locations shown 
on the CAD model. 






























of the embedded triangle overlayed on the estimated stiffness 
map of the ex vivo porcine liver. The actual registration 
parameters, the estimated registration parameters and the root 
mean square (RMS) error between the estimated location of 
all the probed points and their true positions are shown in 
Table |n| for all Examples. 
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Fig. 9. (a) Registration results for Example 3, (b) Estimated stiffness map 

for Example 3 with 119 probed points. The algorithm starts with an initial 
grid of 19 points: 4 corners and 15 uniformly spaced points. 



Fig. 10. (a) Registration results for Example 4, (b) Estimated stiffness map 
for Example 4 with 79 probed points. The initial set for Example 4 consists 
of 7 samples that define the boundary of ROI and 12 uniformly distributed 
samples inside the ROI. 


TABLE II 

Registration Results 


Example x(mm) y(mm) z(mm) ^a;(deg) 6y(dQg) ^^(deg) RMS(mm) 


Actual 

1-3 

5 

10 

-15 

11.46 

-11.46 

5.73 

- 

Guess 

1-3 

0 

0 

0 

0 

0 

0 

- 

CMU 

1 

4.51 

9.56 

-15.03 

12.57 

-11.42 

5.53 

0.91 

CMU 

2 

4.31 

9.55 

-14.97 

12.53 

-11.42 

5.53 

1.04 

CMU 

3 

4.56 

9.24 

-14.97 

12.58 

-11.27 

5.69 

1.14 

Actual 4 

5 

7 

-13 

11.45 

-5.72 

8.59 

- 

Guess 

4 

0 

0 

0 

0 

0 

0 

- 

CMU 

4 

5.78 

6.4 

-13.04 

11.84 

-5.50 

8.66 

0.74 


IV. CONCLUSIONS 

This work introduced a probabilistic estimation of the stiff¬ 
ness distribution of the organ using Gaussian processes and 
Bayesian optimization to direct the probing for maximum 
information gain. We believe fusing intraoperative data into 
the preoporative model is important to alleviate the limited 
situational awareness in MIS. The performance of the method 
was demonstrated by a number of examples and the results 
show that information-guided probing can avoid probing the 
entire organ and successfully reveal the stiff regions while 
registering the tool to the a priori geometric model of the 
organ. 


There are several directions for future work. We used a 
simple experimental setup and an unconstrained environ¬ 
ment to evaluate our method to avoid additional sources 
of error such as workspace limitation and deflection of the 
robot. In our future work, we will demonstrate the proposed 
method in real-time using a continuum robot [20] and the 
da Vinci Research Kit. Another extension of this work is 
the intraoperative reconstruction of the organ surface and 
stiffness features as the organ goes through changes during 
the surgical procedure. We envision that the information- 
guided probing will enable generation of an updated model of 
the visible anatomy and reduce the time it takes to reconstruct 
the intraoperative scene. 


APPENDIX 


A. Complementary Model Update for CARE 

The CMU method is briefly described here for com¬ 
pleteness. The various steps involved in the CMU can be 
described as follows: 


1) Collection: In the collection step, pairs of force-position 
measurements which satisfy the following conditions 
are grouped together in the same set: 

i) The force magnitudes are different. 

ii) The position measurements are spatially close by. 
Hi) The normal directions are the same or similar up 

to a threshold. 


In the experiments described in Section III-A we as¬ 
sume that the surface is smooth and frictionless, thus the 
applied force is along the surface normal and increases 
with depth. The three conditions stated above imply that 
distinct measurements that lie on the surface normal 
experience different force, and form a compatible set. 
The magnitude of the sensed normal forces are denoted 
by Fk G R, and G are the coordinates of the 

sensed points, where [ denotes that the entity is 

in the tool’s reference frame. Given the measurements 
(p^,F/c), k = 1,2,...,/ obtained so far, we collect 
compatible sets, = 1,2, ...,n, where n is 


the total number of distinct sets obtained and j G 
Ni where Ni is the set that contains the indices of 
measurements belonging to the set. 

2) Stiffness estimation: Eor each set i that has at least one 
pair of force-position measurements, we estimate the 
local stiffness q, assuming a linear stiffness model. 


Ci = L (||(p«)i - {p^)4,{{Fp)i - iF^)i)) (8) 

where {ffy) G Ni, (3 ^ j. In Eq. L(depth,force) is 
the function that returns the slope of the best line that 
fits the variation of force with deformation depth. 

3) Correspondence: The sensed points {p^)i are trans¬ 
formed to the CAD frame using the best registration es¬ 
timate, T. Then, we find (pf, nf) = M ^T(p^)^, 0^, 

where = M(a^,0) is the rule that finds the 
closest point 6^ G 0 to and the corresponding 
normal , where [ denotes entities represented in 
the preoperative CAD model’s reference frame. 













4) Minimization: The following objective function is min¬ 
imized using Arun’s method [21] 


argmm 


i=l 


pf - - T{p?), 


. (9) 


5) Upon obtaining T, we loop between the Correspon¬ 
dence and Minimization step until convergence or up 
to a fixed number of iterations. 

The minimization step can return a local minima when 
Arun’s method [21] or update step of a Kalman filter [22] 
is used. Therefore, we seed the algorithm with multiple 
initial guesses to overcome the problem of local minima. The 
detailed derivation of the algorithm is described in [23]. 
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